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Abstract 



A computer method for objective computation of adiabatic three- 
dimensional trajectories in the free atmosphere is presented. This 
technique utilizes a digital adaptation of the graphical iterative method 
introduced by Danielsen (1961). 

The program is applied to a test case and is found to give generally 
satisfactory results. Some difficulties are encountered in cases of very 
low wind speeds and also in regions of strong anticyclonic wind shears . 

A method for inclusion of diabatic effects is outlined in detail. 

This approach may allow investigation of more complicated flow problems 
than have been previously possible. 



1 . Introduction 



The determination of three dimensional trajectories by air parcels 
in the free atmosphere has been of interest to meteorologists for many 
years. Knowledge of parcel trajectories is valuable in analysis of trace 
substance transports, in the diagnosis of large scale systems, and also 
can be of importance in weather forecasting. 

Although the large scale motion of the atmosphere is overwhelmingly 
horizontal, it has long been recognized that neglection of the vertical 
velocity component can produce appreciable trajectory errors in as little 
as 12 hours (e.g., see Petterssen, 1956). However, since the large 
scale vertical motion must be inferred rather than measured directly, 
progress has been very slow in developing satisfactory trajectory techni- 
ques . 

Because diabatic effects in non-precipitating regions are generally 
accepted to be quite small in the free atmosphere, the concept of evalua- 
ting trajectories on isentropic surfaces has long been an appealing one 
(Panofsky, 1947; Saucier, 1955; Petterssen, 1956). Many different 
techniques have been used previously -- the most straightforward being 
a direct graphical calculation from isotach and streamline analyses 
(Palmen and Newton, 1951; Danielsen, 1959; Reiter and Danielsen, 

1960; Staley, 1960, 1962; Reiter, 1963). 

In any trajectory scheme the accuracy of the computed displacement 
tends to degenerate rather seriously with time. This is due to the 
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uncertainty inherent in basic data which is only available in discrete 
form in both space and time. Calculations of trajectories by the above 
method directly from streamline and isotach analyses are particularly 
difficult due to the extreme sensitivity of the calculation to errors in 
the streamline analyses. For example, a wind direction error of just 
10° can produce a normal displacement from the true trajectory of about 
2° latitude in just 12 hours under average wind conditions. The effects 
of errors of this type and of radiosonde spacing errors were pointed out 
by Djuric (1961) . 

In an attempt to minimize the uncertainty due to streamline errors, 
Danielsen (1961) developed a more objective technique based upon a 
path integration of the total energy equation expressed in isentropic 
coordinates. This basic graphical technique has subsequently been used 
with good results by Danielsen (1964), Danielsen, Bergman, and 
Paulson (1962), Houghton (1962), Mahlman (1965), and by Reiter and 
Mahlman (1965a, b, c) . This approach has been used in tracing move- 
ment of radioactive debris , with computed results in agreement with 
observed transports. 

2 . Methods of computing three-dimensional trajectories 

As noted in the previous section the energy method developed by 
Danielsen (1961) has proven to be an effective technique for obtaining 
three dimensional trajectories in the free atmosphere. A basic diffi- 
culty in employing this system is the time of computation involved. 
Initially, isentropic surfaces must be generated containing analyzed 
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wind and Montgomery stream function (c T + gZ) information. An excel- 

P 

Lent treatment of the problems involved in obtaining Montgomery stream 
functions on isentropic surfaces is given by Danielsen (1959). Next, in 
the evaluation of the trajectories, a graphical iterative scheme is 
employed. This becomes very time consuming and cumbersome when a 
large number of trajectories must be determined. 

One way to avoid this difficulty without loss of computational 
accuracy is to perform the calculations with a digital computer rather 
than by graphical methods. This, of course, is advantageous in that 
any such technique is completely objective, once it is successfully 
programmed. Several previous investigators have developed computer 
techniques for evaluating three dimensional trajectories. 

Nagle and Clark (1966) employed a kinematic scheme on isentropic 
surfaces using the sum of the geostrophic wind and a divergent wind 
inferred from the continuity equation expressed in isentropic coordinates. 
This procedure may not yield the most satisfactory results because ihe 
divergent part of the wind field inferred by this method may be far from 
reliable. Also, Krishnamurti (1966a) showed that appreciable ageostrophic 
wind components are mostly non-divergent. This difficulty could possibly 
be reduced by using the balance equation to deduce the non-divergent 
stream function. 

Paegle (1966) computes three-dimensional trajectories kinematically 
in a pressure coordinate system using a vertical motion (co) in the 



3 



vertical and the sum of a non-divergent and divergent wind in the hori- 



zontal. The non-divergent stream function (ip) is obtained from the non- 
linear form of the balance equation. The velocity potential (\) is ob- 
tained from the vertical motion field. This vertical motion field is derived 
by a succession of iterations of the complete non-linear omega equation 
(Krishnamurti , 1966b). This technique contains the possibility of giving 
excellent results provided that the derived three dimensional winds are 
not overly sensitive to errors in the initial height field and further, are 
not excessively smoothed during the complex numerical procedures. In 
fact, part of the objective of this study is to develop an alternative 
trajectory method to this one in order to provide more rigorous tests 
for each approach . This is particularly important, since both schemes 
are very complex numerically, and require considerable preparation of 
input data . 

In the usual kinematic trajectory scheme one integrates the equations 



d 0 = — dt 
a 



( 1 ) 



dX = 



u 



dt 



a cos 0 

where 0 is the latitude, X the longitude, u and v the horizontal 

scaler velocity components, and a is the radius of the earth. Because the 

u and v velocity components are either explicitly or implicitly derived 

from the analyses of streamlines and isotachs, errors in streamlines 

can produce appreciable errors in the final trajectory positions. In the 

free atmosphere the wind shears normal to the flow are generally much 
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larger than those along the flow. As a result of this, streamline errors 
producing computed trajectory positions normal to the true position are 
especially serious. 

A way to reduce the problems associated with this difficulty would 
be to devise an objective method which "constrains" the parcel from 
deviating so far laterally from its true position. The graphical method 
proposed by Danielsen (1961) is specifically designed to provide this 
constraint by using the additional information given by the distribution 
of thermodynamic variables. 

Because this graphical method has produced highly reliable results, 
a proper computer simulation of this technique should also prove to be 
reliable. The rest of this section will discuss the concept behind this 
method and demonstrate how it is done on a computer rather than graphi- 
cally. 

Since this method is so fundamental to the subsequent development 
of the machine method, the basic concepts of this system will be pre- 
sented here . 

The horizontal equation of motion in isentropic coordinates is 
dV , _ _ 

— =- fk xV 2 -v e M + F (2) 

where is the horizontal wind vector, f the Coriolis parameter, k 
a unit vector in the vertical, Vq the two dimensional gradient operator 
on a surface of constant potential temperature 6, M the Montgomery 



5 



stream function, and F is the friction force. 

Taking the dot product of the horizontal wind vector into Eq . 
gives 



_d 

dt 





(2) 

(3) 



the mechanical energy equation expressed in 0 coordinates. Performing an 

dM 9M 

Eulerian expansion on the first term and noting that - V- • 7 M =- - — — + t— , 

2 e dt st e 



V • V V ' V 

d / _2 2\ d0 A/_2 2 \ _ dM SM r* . ? 

3V^V+dt zey 2 J''dt„ + MT +V 2 F - 

0 0 0 



(4) 



The substantial derivatives with 0 subscripts now represent the 
substantial derivative computed without the vertical advection term. 

If Eq . (4) is now integrated over the time interval /it, one obtains 




(5). 



-f f Te (^)de + j'v-Fdt, 

e i l i 

where the i and f subscripts represent the initial and final values, re- 
2 -> -» 2 2 

spectively^nd V = v 2 * = u +v . The left hand side of Eq . (5) 

represents the change of the Montgomery stream function over the time 
interval of the projection of the true trajectory onto the 0 surface. 

The first term on the right hand side of Eq . (5) (term 2) is the change 
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of kinetic energy of the projected trajectory. Term (3) represents the 
integrated local M tendency along the trajectory path. Term (4) is the 
change of kinetic energy on the original 0 surface due to non-adiabatic 
motions. Finally, term (5) is the work done on the mass of air by friction. 

Eq . (5) is now in a form such that it can be used to produce a better 
approximation to the final trajectory point. In general, this equation 
must be solved iteratively since the final value of M determined from 
it is a function of the true trajectory path. Eq . (5) thus serves to provide 
a constraint on the deviations of a computed trajectory normal to its true 
position resulting from streamline errors. Since the Montgomery stream 
function on a 0 surface is exactly analogous to the geopotential height 
on a pressure surface, a restraint on the permissible values of 
provides a satisfactory method for reducing the error of a computed tra- 
jectory . 

In graphically computed isentropic trajectories the usual procedure 
is to determine the first approximation to the true trajectory by evalua- 
ting Eq . (1) in the integrated form 

S = V At , (6) 

where S is the displacement and V is the average wind speed along 
the trajectory over the interval At . The information from this initial 
approximation is then used to compute a value of from Eq . (5) 
with terms 4 and 5 omitted. If the computed value of from Eq . (5) 
differs from the value of measured from the initial trajectory, the 
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initial trajectory is then displaced laterally until the final point on the 



chart satisfies the computed value of M^. At this point Eq . (6) is re- 
evaluated along the new path to correct for possible changes in wind 
speeds along this new path. This new trajectory is then used to re- 
compute the terms in Eq . (o) . If the computed and actual values agree 

7 

within predetermined limits (usually ~ .1 x 10 erg/g ), the trajectory 
is completed. If not, the process is repeated until the desired accuracy 
is reached . 

Even though, in general, air movement in the free atmosphere is 
neither frictionless or adiabatic, it is commonly accepted that these 
assumptions are approximately valid over short periods of time. According 
to Kung (1967) a reasonable value of energy dissipation in the free atmos- 
phere is 2.3 erg g ^ sec * or .01 x 10^ erg g * 12 hr If this value 

7 -1 -1 

is compared with the typical values of .5 x 10 erg g 12hr for terms 
1,2, and 3 in Eq . (5), its neglection is clearly justified. However, in 
the surface layer its effect becomes more important. The effect of 
neglecting diabatic motions is at times more serious. Inclusion of these 
effects will be discussed in detail in section 4. 

3. Machine adaptation of the energy integral method 

The previous section has shown how the energy equation expressed 
in 6 coordinates can be used to obtain more accurate three-dimensional 
trajectories in the atmosphere. Because of its utility this system has 
been adapted for solution on the computer by digital rather than graphical 
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methods. Due to the requirement of an iterative solution to this problem, 



there exists an almost unlimited number of approaches which could prove 
successful. The method employed here will be described in some detail 
to illustrate one solution to the problem . 

Although the longer range goal is to apply the technique to the 
generalized problem including diabatic flows (see section 4) the initial 
programming is restricted to the case of isentropic flow. This produces 
an appreciable simplification in the original programming logic. 

All calculations are performed on a Cartesian grid of arbitrary size 
and spacing. The spacing between grid points is fixed in degrees lati- 
tude and longitude. Because of this, the actual spacing in the east- 
west direction is a function of latitude. Basic constants read in are 
the latitude- longitude spacing between grid points, the southern lati- 
tude and western longitude boundaries, number of grid points in the x 
and y directions, number of levels in the vertical, and number of time 
steps to be calculated. The initial variables read in are wind direction, 
wind speed, pressure, and Montgomery stream function at each point in 
the 4-dimensional (x, y, 0, t) grid. 

The first computational step is to convert the velocity field to u 
and v components. Also, the variation of the grid point spacing along 
the x-axis and of the Coriolis parameter with increasing latitude is 
determined. Next, the u and v fields are linearly interpolated to give 
values at times between the 12 hour observation times used for the 
input data. A simple difference is taken to obtain the local 12 hour 
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Montgomery stream function tendency. At present this 12 hour difference 
is assumed to be representative of the tendency over the entire 12 hours. 
In future programs higher order interpolation schemes will be employed. 
This will be at considerable difficulty due to the additional data tabu- 
lations necessary and increased machine storage requirements. 

The first approximation to the final trajectory is obtained by a 
rather straightforward kinematic approach. Beginning at a grid point 
the parcel is displaced for one hour using the initial values of u and 
v to determine the length and direction of the displacement. At this 
point the time interpolated u and v values are spatially interpolated 
from the nearest four grid points to the parcel to obtain u and v wind 
components for the parcel at its new position. These new u and v 
values are then used to extrapolate the parcel for another hour, at which 
time new u and v components are again found. This process is 
repeated until the 12 hour time is reached. 

At each point along the trajectory a test is made to determine 
whether the parcel has left the original grid. If the parcel is found 
to be outside this data field, the trajectory is terminated and a new tra- 
jectory is begun at the next starting grid point. If the parcel remains 
within the grid, preparations are begun for the iterative phase of the 
computation using the adiabatic, frictionless form of Eq . (5). Values of 
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kinetic energy and Montgomery stream function are calculated at the 
initial and final points of the 12 hour trajectory. Also, the finite dif- 
ference form of the local Montgomery stream function tendency is in- 
terpolated at each point and then summed to get a mean tendency along 
the 12 hour trajectory . 

It might be noted that this initial approximation to the trajectory 
can provide a valid result as long as the original wind fields have been 
adequately specified and the time interpolations are satisfactory. As 
a separate problem it is of interest to note the differences between 
these original, kinematically determined trajectories and the final ones 
which also satisfy the energy equation. 

To begin the iteration the difference is obtained between the value 

of M^ computed from Eq . (5) with terms (4) and (5) omitted and M 

obtained from the actual trajectory (6M) . If the absolute value of this 

7 - 1 

difference is less than .05 x 10 erg g , the initial trajectory is 
assumed to be correct and the calculation proceeds to a new trajectory. 

If this is not the case, the trajectory is then corrected in the direction of 
reducing this difference. This is accomplished by noting that in 0 
coordinates , 



V 

g 



1 

f Sn , 



(7) 



where V is the total geostrophic wind and n is the coordinate 

g 

directed normal to the M lines. A "correction distance" can be deter- 



mined from (7) in the finite difference form 



( 8 ) 



6n = 



- 6M 
fV 



Note that V , the actual wind, is used rather than the geostrophic 
s 

wind here. This value of 6n obtained from Eq . (8) must then be trans- 
lated into X and 0 displacements so that the exact direction of the 
iteration can be determined. Since the correction must be made normal 
to the trajectory path, this is accomplished by the expressions 



6X = 



5n 



a cos 0 V 



60 = 



6n u 



V ' 
s 



(9) 



where u, v, and V g are evaluated at the final point of the first approxi- 
mation to the trajectory. 

One may note that this technique does not completely eliminate the 
difference between the two values of unless the flow is geostrophic. 
However, the flow may deviate far from the geostrophic and the method 
will still move the trajectory in the proper direction for convergence. 

If the angle between streamlines and Montgomery stream function con- 
tours is greater than 90°, the method would fail. Since this condition 
is extremely rare in the free atmosphere, the method proves to be 
generally effective. 

Once 6X and 60 have been determined for the final point, each 

point along the trajectory is displaced a distance proportional to 6X 

and 60 depending upon its time of occurrence along the 12 hour 

1 9 

trajectory (eg., 6X X hr = 6X, 6X g hf = TT 

At this point a new approximation to the trajectory has been completed. 
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However, this new approximation to the trajectory may be in a somewhat 
different wind speed regime, thus necessitating an additional correction. 
This is accomplished by writing 



5s = S . - S , 

wind geom 



( 10 ) 



where 6s is the distance correction (to be made parallel to the trajectory) 
Swind t ^ ie total displacement along the new trajectory obtained by 
summing 1 hour displacements using the total wind at each point. S 

geom 

is the distance obtained by summing the actual geometric distances be- 
tween the 1 hour trajectory points. Thus 6s provides a correction for the 
actual wind speed in the vicinity of the new trajectory. Note that this is 
only a correction for the speed and does not correct for direction. This 
value of 6s is then converted to X and 0 displacements according to 

6X=--^— £ 60 = — ^ (11) 

a cos 0 V a V 

s s 

At this point the trajectory has been corrected for the Montgomery 

stream function change plus an additional correction for distance. The 

difference between computed from the adiabatic, frictionless form of 

Eq. (5) along the new trajectory and the value of actually measured 

from the trajectory end point is then recomputed. If the absolute value 

7 -1 

of this difference is now less than .05 x 10 erg g , the trajectory is 
complete. If not, the iterative cycle is repeated until convergence is 
reached . 

Normally, the iterative correction for a 12 hour trajectory is not 
particularly large (usually less than 2° latitude). However, in certain 
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situations when appreciable accelerations are present, or the large 

scale systems are moving rapidly, the corrections can be substantial. 

7 - 1 

Fig. 1 is an analysis of Montgomery stream function (10 erg g ) 
for a trajectory test case at 6 = 300K, beginning 22 November 1962 , 

1200 GMT. Pressures at 6= 300K are given in mb. 

A sample of some of the adiabatic trajectories obtained from the 
machine method outlined above are given in Fig. 2 for the 12 hour period 
beginning on 22 November 1962 , 1200 GMT. The dashed lines are the 
kinematic trajectories used as a first approximation to the final tra- 
jectories. The solid lines represent the 12 hour trajectories obtained by 
the energy integral method. A single solid line emanating from a point 
means that the kinematic trajectory also satisfied the energy integral 
condition. The two types of trajectories tend to be similar — the 
iterated trajectory normally characterized by a small lateral correction 
to the kinematic one. 

A marked exception to this general characteristic may be seen in the 
lower left trajectory in Fig. 2. In this case the iterated trajectory is 
displaced almost perpendicularly from the kinematic trajectory. The 
cause of this difficulty may be seen in Fig. 1, which shows that this is 
an area of very weak M gradients. Because the lateral distance cor- 
rection is made from Eq. (8), a region of very low wind speeds may lead 
to an unrealistic correction distance. Further, because the distance 
of correction is relatively large for a given value of 6M, in such cases. 
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FIG. 1. Analysis of Montgomery stream function (10 erg g ) at 0= 300K for 22 November 1962, 
1200 GMT. Pressures at 0= 300K are plotted. Enclosed region is area of trajectory computations. 
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Fig. 2. Selected isentropic trajectories at 6 = 300 K from 22 November 1200 GMT to 23 November 
0000GMT, 1962. Dashed lines represent the kinematic trajectories computed directly from the wind field. 
Solid lines are the energy iterated trajectories. If no dashed line is present, the kinematic trajectory also 
satisfies the integral condition of Eq . 5 . Pressures are given at initial and final trajectory points . 



the correction is particularly sensitive to analysis errors in the M fields. 
Thus in this case, the "corrected" trajectory becomes less accurate than 
the kinematic trajectory. In ihis program, and also in graphical calcu- 
lations, it has been found to be advantageous to simply use the kinematic 
trajectory in cases where the wind speed is less than 10 m sec * . 

Another source of error in the iterative calculations arises when 
there exists a very rapid change of wind direction perpendicular to the 
estimated trajectory path. In such cases the energy iterated trajectory 
direction may differ somewhat from the actual wind direction at the final 
point. This type of error also usually occurs when the wind speed it- 
self is quite low. As above, the kinematic trajectory is then held to be 
the more accurate one. 

A final difficulty in the iterative scheme arises when the wind shear 
is strongly anticyclonic in the region of the trajectory. In this case the 
changes of Montgomery stream function and of kinetic energy are nega- 
tively correlated in the iterative correction of an approximate trajectory 
using Eq . (5). For example, a correction to increase the value of 
will simultaneously cause a marked reduction in the final kinetic energy. 

A result the computation may be left as far from convergence as in the previous 
approximation. Occasionally the iteration may even diverge from the true 
solution. In such cases the kinematic trajectory is again used as the 
better approximation to the true trajectory. 

In regions where the above difficulties are not present, the iterative 
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procedure gives very acceptable results for adiabatic trajectory calcu- 
lations . The inclusion of diabatic effects produces considerable 
additional complication and will be discussed in the subsequent section. 
4 . Incorporation of diabatic effects into trajectory determinations 

Although the air flow in the free atmosphere tends to be roughly 
adiabatic in character, diabatic effects are always present. At times 
these diabatic effects may lead to a true trajectory which is appreciably 
different than that obtained by assuming isentropic conditions. 

At first it might seem to be contradictory to consider diabatic effects 
within the framework of a problem formulated in isentropic coordinates . 
However, even for diabatic trajectories, the 0 coordinate system retains 
its advantage over Z or p coordinates for several reasons. First, the 
adiabatic component of vertical motion is invariably larger than the 
diabatic part. This implies that the relative distance a parcel will 
"depart" from its original isentropic surface is appreciably smaller 
than the departure from a Z or p surface. Second, the diabatic effect 
of release of latent heat of a parcel is strongly dependent upon the 
cooling of an air parcel due to adiabatic ascent. Thus, an accurate 
determination of the three dimensional motion of a parcel is logically 
obtained by utilizing an initial knowledge of the adiabatic part of the 
vertical motion. 

This fact will now be used to develop a method for evaluating the 
diabatic components of the motion as an extension of the isentropic 
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method outlined in previous sections. From the definition of a mean 



value , 



* • 2 dp 

= J 5T 

'i 



t 1 ad 



, r> 2 dp , 

, dt + , dt 

dt, J dt, . 
t^ lr t^ lh 



( 12 ) 



where co= dp/dt, dp/dt^ is the adiabatic part of to, dp/dt^ is the 
contribution to to from long wave cooling, and dp/dt^ is the to due 
to release of latent heat of condensation. 

However, the second and third terms on the right hand side of 
Eq . (12) represent a pressure change due only to a change of potential 
temperature of the parcel. Thus one may write, 

“ At = A e a d + f A 8 lr + H Mh < 13) 

where 3p/S 0 is the measure of the mean static stability along the tra- 
jectory and A0 is the change of potential temperature. Thus, in this 
model the total vertical displacement of an air parcel is due to a combi- 
nation of three effects, only the first of which is available from the 
method described previously. The term in Eq . (13) involving A0j r can 
be extremely complex to evaluate, since it, in general, includes long 
wave radiative effects of water vapor, ozone and carbon dioxide. How- 
ever, the contribution of ozone is very small in the troposphere. Also, 
according to Manabe and Moller (1961) the contribution of carbon 
dioxide is considerably smaller than that of water vapor. 

Making use of the above fact, Danard (1968) has developed a 
technique for computing the long wave cooling due to an arbitrary cloud 
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and moisture distribution. Danard's model has been adapted for appli- 
cation to this problem. Fig. 3 shows diabatic cooling rates (deg. C/12 hr) 
computed from Danard's model and interpolated to the 0 = 300K surface 
for 22 November 1200 GMT. This figure shows that the spatial variation 
of diabatic cooling on a 0 surface is surprisingly large. Cooling rates 
in excess of 2 deg/12 hr are observed in the northeast corner. These 
high values are due to radiation from the top of the cloud cover on the 
east side of the trough observed in Fig. 1 . Note also the low cooling 
rates in the central part of Fig. 3. This area is within and behind the 
base of the trough and has been shown to be a large mass of descending 
dry stratospheric air in a previous paper (Reiter and Mahlman, 1965c) . 
Another reason for the pronounced variability of the cooling rates is due 
to the large slope of the 0 = 300K surface. In this case the pressure 
varies from 920 to 400 mb, thus leading to large variations of height 
and temperature over the map. It may be seen, however, that this 
method contains the possibility of gaining a more adequate determination 
of the diabatic cooling of a parcel due to long wave radiation. The most 
straightforward way to include this in a trajectory calculation would be 
to simply use the average of the cooling rates computed at the initial 
and final points of the adiabatic trajectory. 

The last term in Eq . (13) represents the contribution to the total 

pressure change of a parcel due to release of latent heat of condensation. 

This, of course, is a heating effect and tends to balance the diabatic 

cooling due to long wave radiation. The heating of an air parcel due 

20 
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FIG. 3. Computed diabatic cooling rates for 22 November 1962 , 1200 GMT at 0= 300K. Units are 
expressed in deg C/12 hr. Strong cooling rates at eastern edge are due to cloud top radiation. Low cooling 
rates in central region are associated with the dry descending air in this region. 



to condensation of water vapor may be written with sufficient approxi- 
mation as 



dH _ C p T dQ 
dt ih" 6 dt ih 



Sq s 

Luj-^- 6 (x , y , 0, t) 



(14) 



where — — is the heating due to release of latent heat, L the latent heat 
dt lh 

of vaporization, Bq^/Bp the rate of change of saturation specific 
humidity with respect to pressure and 6(x,y,0,t) is either one or zero 
depending upon whether the point x, y, 0, t is saturated or not. 

To obtain the actual value of heating of a parcel over a 12 hour 
time interval due to latent heat release, Eq . (14) must be integrated 
over At’ , the interval of time which the parcel is saturated. This 

I 

value of At is determined by first determining the initial value of q 

for the parcel and then finding the pressure (p ) at the 0 surface 

sat 

for which the value of q equals the initial value of q. From this one 

s 

may define 



Ap “ p i - P sat (15) 

I 

where A p is the pressure change a parcel must undergo to experience 
saturation, and p^ is the initial pressure. Eq . (15) may then be used 
to form the expression 



At 





(Ap 



ad 



_A E 

+ H A ^r ) 



'N 

j ! 



(16) 



where, if At' 



is less than zero, it is set to be zero for the computation. 
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Here ' < Ap ad + Is A V 

experienced by the parcel 



represents the total 12 hour pressure change 
due to the effects of adiabatic motion and long 



wave radiation. 



One may now integrate Eq . (14) over the time interval At to obtain 

(16a) 
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A „ p f 6 L to s 1 

Ae n. = " J — dt, or 
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Ae ih = 



-1000 X L 
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where x = 0.2857 . The approximation in Eq . (16b) is justified because 

all terms involving deviations from the mean in the expansion of 

w 3q /9p/p are more than one order of magnitude smaller than the term 
s 

retained . 

If the expression for A8^ in • (16b) is substituted into Eq . (13), 
one obtains 
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Note that the o? values obtained from the adiabatic trajectory and the 
long wave cooling effect are used as the initial approximation to the 
true 0) field required in the third term. The subscript on oT^ denotes that 
this is the first approximation to the total pressure change by the parcel. 

Eq . (17) can now be solved iteratively for the final mean 03 by writing 
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where the n subscript indicates the number of iterations. This iteration 
can be continued until I oo , - cc I < c, where e is some acceptable 
error. The mean vertical motion over the 12 hour period has now been 
obtained. Now that the total vertical displacement has been obtained, 
the total diabatic heating of a parcel is given by the expression 
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where A0 determined along the trajectory. An interpolation in 0 is 
utilized to find wind speeds at each new potential temperature along the 
path. Utilizing these new values, a correction is then made for distance 
according to Eq . (10). 

To begin preparation for the lateral correction of the diabatic trajectory, 
A 8 from Eq. (19) is substituted into term (4) of Eq . (5) written as 
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The vertical derivative of the kinetic energy is evaluated along the tra- 
jectory by a centered finite difference of the kinetic energies of the levels 
above and below the parcel at a given point along the trajectory. 

The computed value from Eq . (20) is then substituted into Eq. (5) to 
obtain a new computed value of for the trajectory. A lateral correction 
is then made for the trajectory in the same manner as done previously for 

the adiabatic case. This correction is normally small because the numerical 

7 -1 

value obtained from Eq . (2) is nearly always less than .2 x 10 erg g 



considerably smaller than the usual magnitude of terms (1), (2), or (3) in 
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Eq. (5). 



If this correction is small enough that the computed and observed 
values of agree within the predetermined limits, the trajectory is 
complete and simultaneously satisfies the energy equation and the 
diabatic effects. 

Provided that a correction is necessary in order to satisfy Eq . (5), 
the next step is to provide an additional correction for the change of 
distance along the trajectory due to effects of horizontal wind shear. 

This corrected trajectory is then checked by a recomputation of Eq . (5). 

If the computed and observed values of (projected from the diabatic 

trajectory onto the original 0 surface) do not agree within the predeter- 
mined limits, this process is repeated until the convergence criterion is 
satisfied . 

However, even though Eqs . (5) and (6) are now satisfied for this new 

trajectory, the possibility remains that the computed vertical displacement 

is again in error, since a path correction has been made. A new check is 

accomplished by recycling through the process outlined from Eq . (14) 

through Eq. (18). If oi , is sufficiently close to the previous value 

n+i 

of a^ +1 obtained from the first cycle, the trajectory is complete. If these 
two computed values differ appreciably, the cyclic procedure beginning with 
Eq . (19) is repeated until one of the exit criteria is satisfied. 

At the present time little quantitative information is available as to 
the number of cycles required for the computed trajectory to satisfy the 
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convergence criteria . In a set of graphical computations of diabatic effects 
on trajectories by the author in a previous publication (Reiter et. al. , 

1965), only one to three iterations were normally required. In that com- 
putation, conservation of equivalent potential temperature was imposed 
as an additional criterion to be satisfied along the trajectory. 

Although there are no fundamental difficulties in programming the 
inclusion of diabatic effects for machine computation, in actual practice 
this will prove to be very difficult and time consuming. Further, it should 
be noted that this method will be incapable of determining trajectories in 
areas dominated by sub-synoptic scale motions. 

5 . Summary 

A working computer method has been outlined for objective computation 
of isentropic trajectories in the free atmosphere using a modification of 
the scheme originally derived by Danielsen (1961). This method has been 
found to give highly acceptable results in most cases. The method was 
found to give poor results in either areas of very low wind speeds or 
regions characterized by strong anticyclonic wind shears . 

A general technique for expanding this method for inclusion of diabatic 
flows has been presented in some detail. This should enable investigation 
of more general three-dimensional flow problems than have been capable 
of solution up to the present time. The required expansions of the 
isentropic trajectory program to include diabatic effects are currently in 



progress . 
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